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Abstract 

Process zone at the tip of a propagating crack engendered by the stress-induced 
local phase transition of the second order is studied theoretically. We show that the 
zone can only exist within a certain domain of the phase diagram at one side of the 
phase transition line depending upon the sign of the striction constant. We obtain the 
boundary of this domain and establish its dependence upon the crack velocity. We 
show the existence of a critical crack velocity above which the zone cannot exist. We 
report the overcritical solution for the order parameter describing the incipient process 
zone, while far from the bifurcation point we solve the problem numerically. 


1 Introduction 

It is important to study the brittle fracture both from the point of view of numerous ap¬ 
plications and in order to further propel our knowledge of intrinsic properties of solids. A 
set of experimental facts presently collected on fracture is very impressive. However, the 
understanding of physical aspects of the its mechanisms is still far from completeness. 

Since decades a common opinion worked out that the crack behavior is determined by 
the process zone, a nano- to mesosized domain in the vicinity of its tip. For the long time, 
however, the process zone was beyond the reach of experiments. Recently new experimental 
techniques (described in more details below) emerged that enable one to study structure 
and properties of the solid in the close vicinity of the crack tip. These as well as some 
more traditional experimental techniques have shown that fracture is often followed by a 
rearrangement of the solid structure in the close vicinity of the crack tip. This challenges 
one to hnd a theoretical approach adequately describing such a phenomenon. 

Stress-induced local phase transitions (LPTs) at crack tips have been reported in litera¬ 
ture since long time. They have been observed in different classes of materials. 

Stress induced austenite-martensite LPTs in metals have been studied since 60th [5] and 
are under a keen attention up to now. Martensite LPTs have been observed in iron [6] and 
steels 0. i- Active studies of the martensite LPT in Ni-Ti (Nitinol) alloys widely-used in 
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application are being actively carried out 0. cni. El. LPTs are also exhibited by other 
shape-memory alloys, such as Cu-Al-Ni [12], [13], Cr-Ni [H], Ni-Al [15] and Ti-Al-Nb [T6] . 

Crack tip LPT from the bulk bcc phase into the nanoscale fee phase zone has been very 
recently observed in molibdenum na, the latter phase being nonexistent under pressure, 

p > 0. 

Zr 02 based ceramics have received a large attention, since it has been observed that 
the tetragonal-monoclinic phase transition at the crack tip strongly improves their fracture 
toughness [TS], [20] [21] • Ferro- and antiferroelectric ceramics have also been reported to 
exhibit LPT under fracture [22], [I]- 

The fracture toughness improvement due to the superconducting LPT in YBCO and 
BSCCO at the crack tip has been reported in the paper [2^ . 

A structural rearrangement within the crack tip zone in sapphire manifested in formation 
of metastable Al-O-Al clusters at its fracture surface has been recently reported [23]. 

The crack tip stress-induced structural LPTs have been also observed in polymers and 
epoxies [20] • Resins are known to exhibit crystallization at the crack tip strongly affecting 
the fracture process [26] , 

As for general trends in the evolution of the methods of the LPT observation one can note 
a transition from indirect methods of analysis to those making it possible to obtain a direct 
structure of the transformed zone with a high spatial resolution, combining few techniques 
in one study often appropriate for propagating cracks. One should first of all mention the 
method of high-angle annular dark-held scanning TEM approach allowing for the direct 
imaging of atomic locations. It is this method that has been recently used to detect the 
crack tip LPT in molibdenum. To exclude any interpretation ambiguity the results have been 
further combined with the electron nano diffract ion patterns El- Further, the combination 
of micromechanical loading with in-situ high-resoluion X-ray microdiffraction m enabled 
the authors to image a complex LPT zone in the polycrystalline Nitinol. The combination 
of the in-situ SEM with the electron backscatter diffraction [8] made it possible to study the 
evolution of the emerging phase in the tip vicinity during the fatigue experiments. The in- 
situ optical digital image correlation technique made it possible to obtain strain fields and the 
phase boundaries at the tips of propagating cracks m Raman mapping revealed the local 
distribution of phases in the vicinity of the crack tip [1] . AFM maps the lateral distribution 
of the surface height. The latter is directly related to spontaneous phase transition strain, 
thus, enabling one to distinguish phases and determine the phase boundary 0. E|. 

On the theoretical side three approaches can be pointed out. First, atomistic mechanisms 
of LPTs have been revealed by computer simulations and density functional theory-like 
calculations for several solids, such as iron silicon [20], [31], [32], tantalum [23], 

zirconium [23], UO 2 [23], molibdenum El. Nitinol [23]. This became possible as the result of 
the computation power development, and advancement of the molecular dynamics approach, 
as well as in implementation of hybrid approaches combining the molecular dynamics with 
quantum mechanics [22]- The simulations revealed a strong dependence of the LPT formation 
upon (i) loading mode, (ii) crack plane and direction and (iii) sample geometry [27[, [28] . 
They further elucidated atomistic mechanisms leading to the LPT development [2Ej, [32] . 

Second, a number of researches exploited a mechanical approach treating the LPT zone 
as the one only differing from the rest of the solid by its (i) elastic properties and (ii) 
spontaneous strain. Antolovich [12] was first to propose a mechanism of the transformation 
toughness of a quasi-static crack. References to further papers of this kind one hnds in the 
reviews [19], IZI El as well as in the book [46] . 

The above approach ignores the fact that the local phase transition is related to one 
or several internal degrees of freedom of the solid obeying their intrinsic constitutive laws. 
A third stream of works just focused on the role of the internal degrees of freedom in the 
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formation of the LPT at extended defects. It has been pioneered by the paper of Nabutovsky 
and Shapiro describing the dislocation-induced LPT m- In the papers of Korzhenevskii the 
effect of the LPT on the behavior of dislocations and, thus, on plastic properties of the solid 
has been established |1S] , formation of the LPT at a moving dislocation has been described 
in the paper |19]. LPT at wide domain walls has been described in [50], while those at 
narrow twin boundaries has been addressed in the papers [5T] . This approach has also been 
developed in application to brittle fracture. Formation of the LPT at the crack tip has been 
analytically described in the papers |S2] and [53], and numerically in [51] and [55] . 

It is generally accepted that mechanisms of the brittle solid resistance to the crack prop¬ 
agation are attributed to its process zone (PZ). The PZ notion only has sense, if such a 
zone can be clearly distinguished from the bulk of the solid. The latter can only be done, 
if at least one its physical property exhibits a perceptible variation across the PZ bound¬ 
ary. It may either be an abrupt quantitative variation such as the elastic nonlinearity [H] 
or hyperelasticity [12], or any qualitative variation. In the latter case the PZ differs from 
the bulk by e.g. its chemical composition or crystal structure. Such properties determining 
qualitative differences are always controlled by internal solid degrees of freedom, rj, as e.g., 
concentration of reaction species in the former and phonons, magnons, electronic degrees of 
freedom, etc. in the latter case. 

We focus on the case of the PZ qualitatively different form the solid bulk by its structure 
described. Condensation of any such structural degree of freedom in the bulk of the solid 
{rj = const 7 ^ 0) corresponds to a bulk phase transition (PT): structural, magnetic, electronic 
respectively. The degree of freedom related to the condensate, rj, is referred to as the ’’order 
parameter”. 

High values of stress at the crack tip: a ~ may trigger formation of a transforma¬ 
tional PZ, where r] = ri{r) ^ 0 within the zone, while vanishing outside. This situation may, 
thus, be regarded as a LPT. Here a is the stress tensor and r is the radius-vector counted 
off from the crack tip [39] . 

The present paper reports the case of the second order LPT admitting a fully analytical 
treatment. The latter is of a paramount importance, since it gives hints of what can be 
expected in more complex cases that cannot be treated analytically. 

We show here that in a solid undergoing a second order phase transition a LPT zone 
with the size 

forms at the tip of the motionless as well as propagating crack either above, or below the line 
of the phase transition on the phase diagram. We show that in the both cases the difference, 
T* — Tc, between the transformational process zone (TPZ) emerging temperature, T*, and 
the bulk transition temperature (the Curie point), T^, is 


T* — T- ~ i 


1/3 2t t2 


Aga 


( 2 ) 


The difference, AT* = T* — we refer to as the ’’temperature shift”. Here Kj is the stress 
intensity factor applied to the crack, E is the Young’s modulus, k = dT^/dp is the slope 
of the phase transition line in the (p, T) phase diagram, k is the order parameter kinetic 
constant and g is the one dehning the energy of the order parameter inhomogeneity, the 
parameter a is related to the Curie constant, C as a = 2ti/C. The upper sign corresponds 
to fc < 0, while the lower one - to fc > 0. Therefore, if the phase diagram slope is positive 
{k > 0), the LPT only takes place below the phase transition line in the phase diagram 
AT* < 0, while the negative slope {k < 0) results in AT > 0, that is, the zone containing a 
low-temperature phase is embedded into the matrix of a high-temperature phase ([T]). 
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Figure 1: Schematic view of the location of the process zone on the phase diagram depending 
upon the slopes of the transition lines. (A) The bulk phases I and II are separated by the 
phase transition line, Tq{p) (the solid line), with k > 0. The region of the process zone 
existence, Ila, is below the transition line, where the phase I process zone is embedded into 
the matrix of the bulk phase II. In the region Ilb, below the line T*(p) (the dashed line) 
the process zone vanishes. (B) The case of the negative slope, fc < 0, of the bulk phase 
transition line, Tq{p) (solid). The process zone existence region, la, is situated above the 
phase transition line, where the phase II is embedded into the matrix of the phase I. In the 
region Ib above the T*(p) line (dashed) the zone vanishes. 


We, further, demonstrate the existence of the critical velocity: 

V,r^^ia\k\Kjf^ (3) 

K 

such as soon as the crack tip velocity exceeds the critical one, the TPZ vanishes. 

In the discussion we hrst illustrate the results using the examples of few ferroelectrics: 
BaTiOs, PbTiOa and LiNbOa. We, further, estimate typical values of AT*, valid for most 
inorganic solids. This represents our most striking result. Indeed, the typically AT* values lie 
between ~ 10^ and ~ 10^ K. This implies that the TPZ existence region covers a considerable 
part, if not the whole phase diagram, though only at one side of the phase transition line. 
We then discuss experimental difficulties of the TPZ detecting, as well as possible ways of 
its observation. We, finally, put forward a general approach for the concept of the process 
zone. 

The paper is organized as follows. In Section II we formulate an equation describing the 
order parameter dynamics. Lengthily details of its derivation we give in Appendix A. In 
Section III we hnd the analytical solution of the PZ equations. In Section VI we report our 
numerical results. In Section V we make numerical estimates and give comments generalizing 
our hndings. 


2 Process zone dynamics 

Dynamic behavior of the order parameter splits into several universality classes as it is 
suggested by the Hohenberg-Halperin-Ma scheme |S7|. Within this scheme the most simple 
and simultaneously most often met case is referred to as the model A describing the dynamics 
of a nonconserved order parameter, p, corresponding to a phase transition. 

We focus on the case related to crystal structure variation within the process zone, that 
is, the case of a structural LPT. For PTs of such a type the order parameter generally 
represents a multicomponet object, p = {pi,p 2 , ...,Pn), each component being composed of 
combinations of displacements of the unit crystal cell atoms. Its properties are completely 
determined by the corresponding irreducible representation of the crystal symmetry group in 
the high-symmetry phase of the solid [58], [59], [60]. In the high-symmetry phase p =0, while 
in the bulk low-symmetry phases the components pi {i = 1, 2,..) of the order parameter, are 
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either equal to constants, or to zero. They are independent of spatial coordinates, at least 
within a single domain of a given phase. Such an approach introduced by Landau in 1937 
[5S] enables one to classify all phases and to predict possible bulk phase diagrams [SH], (SH]- 
Here we adapt this approach for the description of the local transformational PZ. 

In this paper we only consider the case of a single-component order parameter (n = 1). 
Being simple it already catches the most important properties of the dynamic LPTs. It 
is also important for another reason. One observes that in solids with PTs described by 
multicomponent order parameters, most of their bulk low-symmetry phases are described 
effectively by one independent parameter. For example, the low-symmetry phases in BaTiOs 
are generally described by the three-component polarization vector rj =(r/i, r/ 2 , "hs)- However, 
in the tetragonal phase one hnds rj =( 0 , 0 , r]), in the rhombohedric phase - rj ={ri, rj, rj) and in 
the orthorhombic phase - rj ={ri, r], 0 ): all these phases are, thus, effectively single-component 
ones. This situation is characteristic also for structural PTs in many other materials. This 
makes the case of the single-component order parameter distinguished. 

2.1 Equation of motion for the order parameter 

Equation describing the order parameter dynamics can be obtained with the help of the 
dissipation function: 

'■"i/d)’'” «> 

and the free energy: 

F = j^{T],eik)dn (5) 

where k is the kinetic constant, t is the time, <F = $( 77 , is the free energy density, Sik is 
the strain tensor and hi is the domain. Since we consider a thin plate case here, hi represents 
a plane and dQ = dxdy. That is, we assign both F and D to the unit solid thickness in the 
z direction. 

Here we address the simplest case of a PT describing by a one-component order param¬ 
eter. In the one-component case the only possible transformations under the action of the 
crystal symmetry group are either 77 —)• 77 or 77 —)■ — 77 . Being invariant with respect to the 
crystal symmetry, the free energy should only contain even functions with respect to 77 [5S] . 
It takes the following form: 


$(77, Eik) = $pt(h) + (6) 

Here ^{ri,eik) is the free energy density, the function denoting its part responsible for 

the PT itself: 


d /T-T ',2 I 2 I 4 /.T'l 

^pt = 2 + 2 ^ 

where g > 0, and /Sq > 0 are the constant parameters of the Landau potential ([S]), while 
a = a{T—Tc), where a > 0 is a constant, Tis the temperature and T^. is the Curie temperature 
and V ?7 is the order parameter gradient. The case of the hrst order transition /^o < 0 will be 
analyzed elsewhere. The strain tensor, Eik, is dehned in a usual way: 

_ 1 , duk ^ 

2 dxk dxi 


Here Ui is the displacement vector. 


5 




is the elastic part of the free energy density. For simplicity we consider here the 
elastically-isotropic case with 




where A and /i are Lame constants [66], where 


( 8 ) 


Eu E 

(l-2cr)(l + cr)’ 2(1+cr) 


(9) 


yields their relations to the Young’s modulus and Poisson’s ratio, a. The latter should not 
be confused with the stress tensor a = aik- It should be mentioned that the constants g, a, 
/9, Tc and A represent the material constants of the solid in question together with E and a. 

Finally, the, so-called, striction constant, A, already introduced above is responsible for 
the interaction between the strain and order parameter fields and can be either positive or 
negative. It should be noted that the form of the interaction term Arfea in (jHI) implies that 
the phase transition only gives rise to the spontaneous dilatation. It is only this case that is 
considered in the present paper. Our final results are more convenient to express in terms 
of the slope of the phase transition line, k = dTc/dp^ on the {p,T) phase diagram directly 
related to the striction constant A. Indeed, one can represent the term ~ in (EllTD as 
a (T - Tc + Aa-^Sii) rf, yielding 


aE 

Equation of motion can be built on the basis of (jSl llj) as follows [SB] : 


( 10 ) 


hF _ 

dp dp ’ dSik 

where d is the variation sign. One obtains the following system of equations: 

= gdAp -[a- 2Aeii {r)]p - (dop^ 

\ doikjdxk = 0 


( 11 ) 


( 12 ) 


where A is the Laplace operator and Oik = d^/dsik is the stress tensor: 


^ik d\£jjdik 2fJj£ik Ap dik 


(13) 


Here d^k is the Kronecker symbol. The last term, Ap'^dik, in the expression fllOp describes the 
spontaneous stress generated by the phase transition. For simplicity we omitted the inertial 
term in the second Eq. (na, which is valid, if E « c, where c is the sound speed. 

Equations flT^ [T3|) represent the complete system describing the dynamics of the trans¬ 
formational PZ. 

One can eliminate the elastic degrees of freedom, eij, form the equations of motion flT^ 
as it is described in details in Appendix A. After their elimination one comes to the single 
equation of motion or the order parameter: 

= g^V - [a + {r)]p - /dp^ (14) 


In contrast to the strain tensor Sa (r) met in Eq. flT^ describing both the field of the tip and 
that generated the PZ, the tensor £)°^(r) only describes the strain field of the ’’undressed” 
tip, i.e., the one without the LPT (?] = 0). It is given by the the well-known fracture theory 
expression [39] : 



(1 -h a) (1 - 2a) Kj 

F(27rr)^A 


cos(6'/2) 


(15) 
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where r and tp are the polar coordinates connted off from the crack tip and a is the Poisson’s 
ratio. The parameter (5 is expressed in terms of (mill]) as follows: 




(16) 


Eqnation flTT|) exhanstively describes the order parameter dynamics within the PZ. It shonld 
be mentioned that ffTT|) can be obtained by the variation procednre flTT]) nsing the dissipation 
fnnction (|1]) and the effective free energy: 


FeS — Fq 




dVL 


(17) 


derived in the Appendix A. 

For the sake of completeness let ns also mention that in the eqnilibrinm {dr]/dt = 0), 
homogeneous {Arj = 0) state, away from the tip = 0) one hnds the bulk phase rj = 0 
also referred to as the ’’mother phase” at T > Tc, while at T < Tc the bulk ’’daughter phase”, 
rj = ± (—takes place [58]. 

Equation ffTTl) represents the model A according to the Hohenberg-Halperin-Ma scheme 
EH. Its solution is demonstrated in Section III. 


3 Analytical analysis of the equation of motion for the 
order parameter 


3.1 The automodel regime 

Assuming the crack tip propagating with the velocity V along the Ox axis and passing to 
the comoving frame, x' = x — Vt, y' = y, one hnds the equation of motion flTTl) in the form: 


gA„ + nvA 


a A B 


cos(6') 


T] — I3rf = 0 


(18) 


where 


B 



a |/i)| (1 + a)Ki > 0 


(19) 


and one chooses the sign ”+”, if A > 0 and ” in the opposite case. Further, r' = 
and the Laplace operator is dehned as A = d‘^/x'‘^ + d‘^/y'‘^. From here on we 
only use the comoving frame and, therefore, omit the primes. 

If one describes a transformational PZ T]{r) ^ 0 embedded into the matrix of the bulk 
phase (cK > 0, r; = 0), localized at the crack tip, (x,?/) = 0, while vanishing away from it, 
the boundary condition takes the form: 


ri{oo) = 0 (20) 

Because its physical origin is related to atomic coordinates which must be limited, the order 
parameter is everywhere hnite: I?]! < oo [M] . 

For the description of the order parameter distribution embedded into the matrix of the 
daughter phase: {a < 0; y = {—a/fdyF ^ g), one needs to use the boundary condition: 

h(oo) = {-a/yyF (21) 

Below to study the problem at hand we employ methods of the bifurcation theory. 
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3.2 Bifurcation theory: a brief review 

For the convenience of the reader let ns hrst shortly recite some key resnlts of theory of 
bifnrcations [BB] which we use in the following argumentation. Let us consider a nonlinear 
equation that can be written in the form: 


L{a)7^ = ( 22 ) 

where L{a) is a linear operator depending upon the parameter a, rj = ri{r) is a dependent 
function and N{ri) is a nonlinear operator, such that iV(0) = 0. In general both L{a) and 
N can be differential, integral or integro-differential operators. In the case considered in our 
further study L{a) is a differential operator, while is a polynomial. 

One can see that equation fl2^ has a trivial solution t] = 0. Assume that it is stable at 
some a > 0. The trivial solution of (12^ becomes unstable, as soon as a reaches a* > 0, 
equal to the hrst eigenvalue, a* = ai, of the linearized equation (12^ : 

L(a„)T„(r) = 0 (23) 

Here are the eigenvalues belonging to the discrete spectrum of the equation fl22)) . if any, 
'h„(r) are their corresponding eigenfunctions and n = 1,2, ...are the natural numbers. In 
analogy with quantum mechanics a* = ai and T*(r) = Ti(r) are referred here to as the 
’’ground state” eigenvalue and eigenfunction. 

In the close vicinity of the bifurcation point one can obtain the asymptotically-exact, 
overcritical solution of Eq. fl22]) in the form of a series in terms of two small parameters: the 
amplitude, and the ’’distance” form the bifurcation point, a — a*. The bifurcation theory 
[68] ensures, however, that its main term always takes the form: 


7](r) ^ ^T^r) + 0(^^) (24) 

Here, ^ is the amplitude to be determined from the nonlinear equation fl22]) . This can be 
done in several ways. Making use of fl25]) and substituting the main term of the solution 
(l2T|l one can, for example, represent fl22ll in the form: ^ x [L{a) — L(q;*)]T* = 
Multiplying scalarly its both parts by T* one Ends: 


L{a) — L{a^) 






(25) 


where we use the notation: (/, g) = J f{r)g{r)dQ for the scalar product of two functions, /(r) 
and g{r), in the Hilbert space. (12^ referred to as a ’’branching equation” [SB], represents a 
nonlinear equation with respect to only valid if ^ is small. Its solution yields the amplitude 
thus, giving simultaneously the solution fl2T|l of the bifurcation problem valid in the close 
vicinity of the bifurcation point. Full details, theorems and their proofs one can End in the 
book of Vainberg and Trenogin [68] . 

We would like to stress, that though the above recipe heavily involves the solution of the 
linear equation (125]) . it represents in fact the solution of the nonlinear equation (122 p . 


3.3 The process zone embedded into the mother phase 

3.3.1 The bifurcation condition 

Let us assume a > 0 and look for the point of instability of the trivial solution rj = 0 
describing the homogeneous mother phase. It can be done using Eq. (ITsp within the approach 
formulated in the previous Section. 




The correspondence between fl22|) and ffT8|) is established as follows: 


L{a)r] = gAr] + kV^ - 


a A B 


cos(6') 


g = 0 


( 26 ) 


NM = I3ri^ ( 27 ) 

Let ns now apply the above recipe of the bifnrcation theory [68] to the case at hand. 
First of all one hnds that at fc > 0 (corresponding to the sign ”+”) eqnation has no 
discrete spectrum meaning that the solution r/ = 0 of the equation flT5]l is stable at a > 0. 

In contrast to that, at fc < 0 (the sign Eq. ( |23|) has a discrete spectrum at a > 0 
implying that Eq. f[T8|l exhibits an instability. Let us consider the latter case. 

At large value of a equation flT8|) has the trivial solution g = 0. Below the bifurcation 
point one hnds a non-trivial solution g = g{r) ^ 0 and (Jjfc(r) = o'j-°^(r) + 0{g^), where g is 
small. According to its dehnition (|7|) a is expressed in terms of temperature: a = a{T — Tc). 
Equation fl5^ yields, thus, the temperature, T*, of LPT at the crack tip. The exact solution 
of the equation fl23|l with L{a) given by is given in the next Section. 


3.3.2 Exact results for the eigenvalue and eigenfunctions of Eq. (I26|) 

Solution of plays as we see an outstanding role dehning both the bifurcation point, 
a* = Oi, and the overcritical solution, fl2T|) of the nonlinear equation (|22|). Let us solve it. 
The linear part of equation fllSp takes the form: 


gA'^n + K,V X d'^n/dx 


Oir 


B 


cos(f)/2) 


rl/2 


= 0 ; 


T„(oo) = 0; |T„| < cx) 


Making the substitution: 


(28) 


T„(r) = exp 



X V’n(r) 


(29) 


one proceeds to the equation in terms of '0n(r): 


gA'tpn - 


O^n + 


4c/ 


-B 


cos(6*/2) 

y/r 


Ipn 


0 


(30) 


The parameter 


Ri = 


(g\^A 

9 

[bJ 2 

a (1 -|- (j)Ki 


(31) 


hxes the characteristic size of the distribution. Passing to dimensionless cylindrical coordi¬ 
nates 6 and p = r/Ri one transforms 0301) into the following equation: 


where (n 


p dp QQ2 



cos(6*/2)\ 

Vp ) 




0 


(32) 


1,2...) represents the eigenvalues of fl52D . They are related to as follows: 


X g4/3 + 4g j 

Eq. fl3^ represents a 2D Schrodinger equation with an anisotropic potential U (p, 9) = 
— cos(0/2)/yCp shown in Fig. HJ 
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Figure 2: Potential U{p,9) = — cos{6/2 )/of the Schrodinger equation fl52]) . Here p^ = 
x/R and py = y/R. (i) indicates the position of the crack tip. 


Let us look for the solution of Eq. fl32|l in the form '^{p,9) = exp{—ez)f{z), where 
z = cos(6*/2). This yields the equation: 

f'{z) - Aezf'{z) - A{e - z)f{z) = 0 

Making use of the transformation / = p{C) exp{z/e) with ( = (z — l/4£^)(£/2)^/^ brings one 
to the equation in terms of (p: 


‘^"(C) - 2C‘^'(C) + = 0; m = 2[1 - {Ae^) (34) 

The latter represents a Hermitian equation with the eigenvalue m only taking non-negative, 
integer, even values: m = 0,2,4... Using fl3T|l one hnds that the condition e: > 0 only fulhlls 
for the ground state solution m = 0 yielding £ = 2“^/^ and p>{C) = const. Now one can return 
to the initial variables and write down the solution of the equation 0321) corresponding to 
m = 0: 


Ai = 


1 


2 ^ 


0.397; Mp) = 


exp + nvp 


cos(6'/2) 


(35) 


The latter is shown in Fig. O (A). Using fl5^ [5S]1 one hnds the ground state eigenvalue, a*: 

1 ^ 2-^2 


a* = 


2^ phs 4g 

yielding hnally the ground eigenvalue fl37|) and the eigenfunction 


(36) 


3.4 The process zone embedded into the mother phase 

3.4.1 Bifurcation point 

The above solution is valid at fc < 0. It gives the ground state eigenvalue: 

o2/3 

0*1 = a(r,i - TJ = ^3 [a \k\ (1 + a)Kif^ - -j- 

and the eigenfunction: 

T*i(r) = exp 


z/i r r 

</AR, 



+ -pr- cos(^^/2) 


(37) 


(38) 
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Figure 3: Dependence of Lf (solid line) and Lf, (dashed line) upon Ui. 


where ui = V/Vd is the dimensionless velocity, and the critical velocity Vd is dehned as 


Vd(Kj) 


24 / 3 ^ 1/3 

TT^/^K 


[a \k\ (1 + 


(39) 


Its physical sense will be discussed below. 

Let us note that Eq. fl371) yields the case T*i — Tc > 0 shown in Fig. [1] (B). 

Assigning the exponent to —1 one hnds the size, Lf, of the order parameter distribution 
in front of the tip: 


Ln = 


22/3 (3 + t^i) + 2 X (2 + 


(1 + z^i)^ 


Ri 


(40) 


At the tip of the motionless crack (z/i = 0) one hnds the order parameter distribution size 
Lfi ^ 9.25i?i decreasing down to Lfi 2.96i?i at vi = 1. The dependence Lfi = Lji(z/i) is 
shown in Fig. [31 


3.4.2 The overcritical solution 

Substituting the obtained eigenfunction fl38|l as well as fi26|l and fiTTI) into the branching 
equation (i25|l one obtains 


I2 (o — 0*1) ^ + ^4/3^^ — 0 

where the factors (n = 2, 4) are the integrals over the whole plane: 


(41) 


In ( ) 




(42) 


depending upon z/i. 

It should be noted that substitution of (1241135]) into the effective free energy, (IT7|) yields 
the LPT free energy: 


FeS — Fq + RI 




(43) 


As expected, its minimization with respect to ^ brings one back to the branching equation 
(ITT]) giving the alternative way to build the overcritical solution. 

The solution of the branching equation has the form: 





P 


1/2 


1 


a > 0*1 
a < ad 


(44) 


11 
















0.8 


■ The simulation results 
— i^"^/;/'^='0.52+0-46v,-0.15vJ 


0.0 0.2 


v,=F/r,2 


Figure 4: The ratio obtained numerically as the function of the dimensionless 

velocity v and its htting by the polynomial flTSl) . 



Figure 5: Spatial distribution of the order parameter, t]{x, y), in the vicinity of the motionless 
(A) and propagating (B) crack. Two blue lines (i) indicate the crack tip position. 


The integrals In{^i) cannot be obtained analytically at ui ^ 0. We calculated their ratio 
numerically by using a standard NIntegrate routine of Mathematica 10.1 EDI employing the 
even-odd subdivision method with the local adaptive strategy. 

The ratio 1^1“^ obtained this way is shown in Fig. H] versus the dimensionless velocity 
Vi. The numerical result can be accurately htted by a simple polynomial: 

/ 

^ ^ 0.52 + 0.46Z/1 - (45) 

hj 

Thus, the asymptotically exact overcritical solution fl2T)l is obtained. 

3.4.3 The critical velocity 

The result fl44D implies that the LPT only takes place at a < a*, while at a > a* the crack 
tip is undressed. Since at a < 0 the whole bulk of the solid transforms into the phase rj 0, 
the domain in which LPT takes place is restricted to the interval: 0 < a < a*. The latter 
result has a very important consequence. Making use of Eq. one hnds that at, E = Wi 
the LPT disappears. In other words the LPT can only exist at the crack tip at E < Vd, 
while at E > Vd it vanishes. This property is fundamental for any LPT, both of the second 
and of the hrst order. It is a direct consequence of the fact that the order parameter has its 
own dynamics exhibiting an intrinsic characteristic time, and that as soon as E > Vd the 
order parameter in front of the crack tip has no time to evolve from y = 0 to rj ^ 

The spatial distribution of the order parameter (|2T)) is shown in Fig. [5l Here (A) shows 
the order parameter in the vicinity of the motionless crack tip, while (B) displays that in 
the case of a propagating crack. This image is obtained with the velocity value V = 0.514i- 
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One finds that the order parameter is highly localized in the vicinity r ^ 10i?i of the tip 
of the motionless crack. In the case of the moving crack tip (Fig. [5]B) the order parameter 
distribution is compressed in front of the tip and stretched in its back with respect to that 
at the motionless one. The length, Lh, of the order parameter distribution behind the tip of 
the propagating crack takes the form: 


Lb 


1 - V/V,i 


(46) 


Behind the crack the length, Lb, diverges, iiV^Vd (Fig. |3]). It should be noted that at 
V —)■ Vci the amplitude, vanishes. 

Let us summarize the results of the present Section. We have shown that in a solid with a 
crack the trivial solution 77 = 0 is stable at high temperatures (T > Tc), but loses its stability 
at the point a = > 0 corresponding to a temperature T*i somewhat higher than that 

of the bulk phase transition: T*i > Tc. This solution describes a region of the phase 77 7 ^ 0 
embedded into the matrix 77 = 0 representing the transformational PZ at the crack tip. In 
terms of temperature its existence is limited to the domain < T < T*i, while in terms 
of velocity to 0 < Id < Wi- The transformational PZ at the tip of the propagating crack is 
deformed with respect to that of the motionless crack: the order parameter distribution is 
compressed in its front, while stretched in its back. 


3.5 The process zone embedded into the daughter phase 


Let us first find a solution for ri{r) in the low-temperature phase {a < 0). We will assume 
in addition that |q;| is large enough, so that where this solution is stable. Observing that at 
r ^ (g/ayR the terms ~ A 77 and ~ drj/dx in flT 8 |l are much smaller than the others and 
neglecting them, one finds the approximate distribution 77 = 770 (r) in the daughter phase: 


770 (r, 6 ') 


1 


a — B 


cos( 6 '/ 2 )\ 


1/2 


(47) 


To be specific, from two solutions of the equation of state we have chosen a positive one. Let 
us now look for the solution perturbation in the form: 77(r) = 770(r) -|- 6ri{r). Substituting it 
into flTK]) one hnds that the term 67]{r) is subjected to the equation with 


L{a)5rj = g/S.5rj + nV 


86 g 

dx 


2 |a| =F I? 


cos( 6 ') 


6g 


(48) 


and 


N{6g) = 3 



By cos( 6 '/ 2 )\ 

7^ ) 


1/2 

577^ -I- ysg^ 


(49) 


In this case at /c < 0 (the sign in 08]) equation fl48ll appears to have no discrete 
spectrum implying that the solution fl47|) is stable. The discrete spectrum indicating the 
instability, however, exists at fc > 0 (the sign ” in 051) . Below we consider this latter case. 

Applying the analysis already described above to the present case one hnds the eigen¬ 
function (l38il . in which instead of Ri one should take a characteristic size R 2 expressed 
as: 


R 2 



. 2/3 7]-F3 

9 

/ “ 25/3 

akKi{l + a) 


2/3 


(50) 


which is smaller than fl5TD by the factor 2 . The bifurcation point has the form: 


0^*2 — 


1 

21/37^2/3^1/3 


[ak{l + a)KiyR + 


8 ^ 


(51) 
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such that a non-trivial solution 6ri{r) = ,^ 2 'h*(r) takes place at 0 > a > a *2 corresponding 
to the phase diagram shown in Fig. [1] (A). Its amplitude, ^ 2 , should be determined form 
the branching equation. Analogously to the previous case the relation fl5T]) gives rise to the 
critical velocity: 

K 2 = [akKj{l + a)f^ (52) 

limiting the existence of the high-temperature PZ in the matrix of the low-temperature 
phase. The latter is by the factor of 2^/^ larger than Vd- 

Since the overcritical solution fl2T)) for Srj is determined by the same eignfunction fl55]l 
as that for the high-temperature case, the distribution, Sri{r) = ^ 2 dz*(r) has the same form 
as that of the order parameter in the high-temperature PZ. The distribution is shown in 
Fig. ([5]). Analogously to the high-temperature case the distribution, Sri{r), at the tip of 
the propagating crack is compressed in its front, stretched backwards and the relation 
holds. The amplitude, ^ 2 , should be determined from the branching equation. 

These further results are easy, but rather cumbersome to derive, and we, therefore, give 
the calculations in Appendix C. 

Let us summarize the findings of this Section. We found that the inhomogeneous solution 
?7o(i’) describing the order parameter distribution in the low-temperature phase is stable at 
a < a *2 <0. At a = q ;*2 < 0 this solution becomes unstable and the solution, ?7o(r) + 
branches decreases off from rjo, the signs of rjQ and 6r] being different. This means that the 
zone with the high-temperature phase, rj = 0, emerges at the crack tip at a = q ;*2 and exists 
within the domain 0 < a < 0*2 and at 0 < IZ < Vc 2 - 


4 Simulation 


4.1 Rescaling 


Based on the bifurcation theory our analytical results are only valid in the close vicinity of 
the bifurcation point. To study the problem far from the bifurcation we simulated solution 
of the equation flTS]) . We report below the simulation of the daughter phase PZ embedded 
into the matrix of the mother phase. 

The equation flT8|l has been rescaled making the variables dimensionless and minimizing 
the control parameters number: x —)■ diXi, y —)• diyi and y{x,y) d 2 u{xi,yi) with the 
scaling factors di = Ri and ^2 = . The rescaled equation ffTSj) takes the form: 


Aim + 2^/V 


du 

dxi 



cos(6'/2) 

^ 1/2 
' 1 


u — 


0 


(53) 


where Ai = d^jdxl + d^/dyf, ri = {xf + dili = dxidyi and the dimensionless 

parameter q is expressed as follows: 


q = 


54/3 


a 


Let us note that in terms of the rescaled parameters (g, u) the analytical expression for 
the bifurcation condition fl37jl takes the form 

l-Z/2 


2 X 2V3 

while the overcritical solution dMHH at q < Qc is expressed as: 


u{ri,e) = 


exp 


-4 [1 — u cos(f^)] + cos 


(64) 


( 66 ) 
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Figure 6: Order parameter distribution obtained by simulations at z/ = 0 and various values 
of qc - q-. (A) 1.5 x 10-^ (B) 6.5 x 10-^ (C) 12 x 10“^ and (D) 16 x 10“^ 


4.2 Results 

We used a pseudo-time stepping approach representing a version of the iteration method. 
Its description along with software technical details and settings are given in the Appendix 
C. The pseudo-time stepping approach converged away from the points: q = 0 and q = qc 
fl5T|) . In practice, we obtained a good convergence at g > 0.1, above the points of the global 
bifurcation. Closer to the point g = 0 we were unable to get an equilibrium solution. In the 
vicinity to the point g = gc the method produced a small regular error discussed in details 
below. Apart from that the method exhibited a good convergence, enabling us to study the 
distribution of the rescaled order parameter, u{xi,yi) in the vicinity of the crack tip. 

Figure [6] shows the distributions of the rescaled order parameter at the tip of a motionless 
crack {u = 0) at the successively increasing values of the difference gc — g > 0 below the 
bifurcation point. As expected the order parameter grows with gc — g, its distribution being 
close to that described by Eq. fl55ll at z/ = 0 (shown in Fig. [5] A). 

To check the proximity of the simulated results to the analytical one we plotted the 
cross-sections of the above solutions by the plane {xi,u) (Fig. [7]). This has been done by 
sampling the solution points from the layer with the thickness 1.5 along the plane yi = 0. 
Because the simulation results exhibit a regular shift of the bifurcation point we plot the 
order parameter versus the distance from the bifurcation, gc — g, rather then g. The solid 
lines in Fig. [7] show the cross-section of the analytical solution fl55|l at the same values of 
gc —g. One can see that close to the bifurcation point the solution obtained by the simulation 
is rather close to the analytical one (Fig. [7]). 

The distribution of the rescaled order parameter, u{xi,yi) at the tip of the propagating 
crack (z/i > 0) at the same value of g is shown in Fig. [HI One can see that the crack motion 
transforms the distribution in several ways. First, with the increasing of the dimensionless 
velocity, z/i = V/Vd, it becomes lower and vanishes as soon as the velocity achieves the 
bifurcation line. Besides, the order parameter distribution is compressed in front of and 
stretched out behind the tip with respect to that at the tip of the motionless crack. 

The latter is in line with the predictions of the analytical solution fl55|l at xi < 0, |/i = 0: 

m(xi, 0) ~ exp \—4r^^^Xi (1 — z/)] ; a:i < 0 (56) 

describing the distribution behind the crack tip. The distribution length diverges as 4^/^/(1 — 
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Figure 7: Cross-section of the order parameter distribution along the plane ?/ = 0 shown at 
X > 0 corresponding to z/ = 0 and various values of q. Dots show the simulation results at 
different values of Qc — q, while the solid line display their httings with the exptession fl55|) . 
The legend for the dots indicate the q values used in the simulation, while that for the solid 
lines yields q values obtained by the htting. 



Figure 8: Order parameter distribution in the vicinity of the tip of the crack propagating in 
the positive direction of x along the line y = 0. All images have been obtained by simupation 
with q = 0.3. (A) shows the motionless case z/ = 0, while (B) and (C) display the case of 
the propagating crack: u = 0.5 (B) and u = 0.58 (C). The image (D) shows the cross-secion 
of the distributions shown in (A-C) by the plane y = 0. 
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yi=viVci 


Figure 9: Dependence of the dimensionless length, Lh/Ri, of the order parameter trail upon 
the dimensionless velocity, ui = V/Vd for several values of the parameter q. For each value 
the data has been taken unltil the bifurcation line has been approached. Dots show the 
simulation result, while the solid line is guiding the eye. 


z/) as 1 / —1 (see also H6|) . The simulation, indeed, looks close to the exponents (Fig. |8]D). To 
check this we htted the backward part of the distributions by the exponents M(a:i < 0, 0) ~ 
exp{xiRi/L). Figure in] displays the dependence of L/Ri upon ui at different values of the 
control parameter q obtained by the simulations. 

It should be noted that the trail length, indeed, increases with the velocity, but does not 
follow the behavior fH6|) predicted by our analytical approach. This is because the analytical 
prediction is only valid very close to the bifurcation point, which is not the case for most of 
the points shown in Fig. IHl 

The plane {ui, q) is divided into two regions. In the region I {q < qc) the transformational 
PZ at the crack tip exists, while it vanishes at the boundary (g = qc) and does not exist 
in the region II (g > qc). To check this we made simulations by hxing the value of ui and 
evaluating the maximal order parameter value, Umax, at different g points. One can see that 
at large g values Umax vanishes, while emerging after a certain threshold. 

The values of these threshold have been extracted from the above data by htting to the 
function: 


Uc 


M) = 


0, g > qcn 

Mo (gen - q<qcn 


(57) 


where uq and qcn are the htting parameters, and g is the variable. The parameter qcn 
represents, therefore, the bifurcation value obtained from the simulations. The subscript ”n” 
stays for ’’numeric”, to distinguish qcn from the analytically obtained bifurcation boundary, 
qc- The dynamic phase diagram obtained this way is shown in Fig. fITT]) . 

One can see that the simulation gives the values of the boundary points that are regularly 
shifted over about 5% upwards with respect to the line obtained analytically. We did not 
succeed to obtain values close to the line gc(z^). This is due to the singularity ~ in the 
equation fl5^ [69] . 


5 Discussion 

5.1 Our findings 

We formulated a natural approach to describe the LPT at the tip of the crack taking into 
account the order parameter, the internal degree of freedom responsible for the phase transi- 
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Figure 10: Dependence of the maximm value of the rescaled order parameter, Umax, upon the 
control parameter q at different values of the dimensionless velocity ui: disks show ui = 0, 
squares - 0.1, hlled diamonds - 0.2, hlled vertex-up triangles - 0.3, vertex-down triangles - 
0.4, open circles - 0.5, open squares - 0.6, open diamonds - 0.7, open vertex-up triangles - 
0.8. The solid lines show their htting to the function fl57|l . 



Figure 11: Dynamic phase diagram divides the plane into the part I in which the 

transformational process zone is present at the crack tip and II without the zone. The dots 
show the results of simulation. They are connected by a dahsed line to guide the eye. The 
solid line shows the theoretical boundary given by the expression flHTjl . 
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tion. The LPT takes place as the response of the solid to the local stress inhomogeneity. This 
snggests that any type of stress concentrator may generate LPT. In addition to the crack 
tips the role of snch stress concentrators play surfaces, inclusions, grain and twin boundaries, 
dislocations and disclinations. It is, further, clear that for any of this type of the concentra¬ 
tors the LPT (i) decreases the total concentrator energy and (ii) introduces its dissipation 
in the process of motion, thus, influencing its static | 18 ], [H] and dynamic 1371 properties. 

Among the LPT generating stress concentrators cracks are the most powerful ones. LPT- 
induced variation of the cracks properties may considerably alter the ability of material to 
resist fracture. This paper focuses on the LPT at the crack tips. The conclusions of the 
present paper are based on our following fundamental hndings. 

First, we have found the exact solution for the point of bifurcation both in the case of 
the motionless and propagating crack. 

Second, we established the existence of the critical crack velocity controlling the disap¬ 
pearance of the LPT zone at high speeds. 

Third, we have found an analytical, asymptotically exact solution of the nonlinear equa¬ 
tion describing the order parameter distribution close to the bifurcation point. 

Fourth, by numeric simulation we obtained the solutions away from the bifurcation point 
which is in line with our analytic solution. 

It should be mentioned here that the present results will be further used for analysis of 
the fast crack propagation published separately. 

5.2 Effect on the crack behavior of the second versus the first 
order transitions 

In this paper we only addressed the phase transitions of the second order. Second order 
phase transitions are generally considered to be soft. This opinion gives rise to the illusion 
that they should only have a negligible effect on the crack dynamics. This point of view 
is erroneous. The softness of the second order phase transitions only means that the solid 
continuously passes from the state rj = 0 to the state // 7 ^ 0 in the transition or bifurcation 
point. This is in contrast to the first order transitions where in the transition point the order 
parameter exhibits a jump. The effect of the LPT on the crack dynamics is, however, not 
related to the order parameter behavior in the transition point. It depends upon absolute 
values the order parameter can achieve within the LPT zone. In the discussion below we will 
argue that the zone is typically very wide. This implies that the PZ order parameter is most 
often in the saturated state. In this respect the second order phase transition exhibits no 
qualitative difference from that of the first order. In the case of structural phase transitions 
the order parameter can be constructed using the atoms displacements from their positions 
in the mother phase. Their saturation values may achieve the values smaller, but comparable 
to the crystal lattice cell dimensions both in the case of the first and second order phase 
transition. 

It should be, further, noted that for one thing the value of the bifurcation point, T*, 
is obtained by the analysis of the linear part, fl25]l . of the nonlinear equation for the order 
parameter. For the other, the answer to the question, whether the PT is of the second, or 
first order depends only upon the structure of the nonlinear terms of the free energy ([7]). 
This suggests that the result ([2]) for T* is valid both for the cases of the first and second order 
transitions. Here we take this assumption as granted. A rigorous proof of this statement 
will be published elsewhere. 
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Table 1: Material constants of BaTiOs, PbTiOa and LiNbOa 



BaTiOs 

PbTiOs 

LiNbOs 

0 

1 

1 

2.1 [71] 

2 EH 

0.7 [76] 

g (xl0“^® cm^) 

1 iHg 

1 IBS] 

0.3 [72] 

Kjc (xl0®ergcm“^/^) 

0.6 ^2.0 [73] 

1.4 [SH] 

1 izai 

K (xl0"^"^s) 

1 H 

1 IHH 

- 

T™(K) 

1898 [83] 

1443 [H] 

1526 [85] 

7’o(K) 

II-III III-IV IV-V 

394 284 200 ^ 

850 [SB] 

1460 [sg 

k (x 10“® K cm^ erg“^) 

-0.8 -0.3 -0.1 [SB] 

-1 mi 

0.018 [ig 


5.3 Estimates 

5.3.1 Material parameters and results for selected materials 

The temperatnre shift valne, AT*, is of a great importance, since it determines the phase 
diagram region where the PZ inflnences the crack behavior. It is interesting, therefore, to 
have numerical values of AT* at least for some materials. As the examples let us consider 
ferroelectrics BaTiOs, PbTiOs and LiNbOs. 

All the three materials belong to a large family of perovskites, LiNbOs exhibiting a 
distorted perovskite structure [82]. 

The high-temperature phases of PbTiOs and BaTiOs exhibit a cubic symmetry. The 
cubic phase in PbTiOs exists above 850 K, while transforming into the tetragonal one at 
lower temperatures m- In contrast to that BaTiOs exhibits the cubic phase (II), tetragonal 
(III), rhomboherdic (IV) and orthorhombic phase (V). The numeration of the phases is given 
according to [HH]- The phase diagram of BaTiOs can be found in the paper [75|. All the 
transitions of BaTiOs and PbTiOs are of the hrst order. 

The Landau theory enables one to describe possible transitions in the both materials as 
generated by the symmetry lowering of the cubic mother phase, the order parameter being a 
3D vector, ( 771 , 772 , hs) associated with the polarization. In the cubic phase 77 = (0, 0, 0). In all 
the low-temperature phases showing up in the phase diagrams the order parameter has only 
one independent component: ( 77 , 0 , 0 ) in the tetragonal, ( 77 , 77 , 77 ) in the rhombohedric, while 
the orthorhombic phase is described by ( 77 , 77 , 0 ) [59] . Effectively, therefore, the problem is 
reduced to a single order parameter enabling one to directly apply the approach developed 
in the present paper. It should be noted that the free energy of BaTiOs and PbTiOs has 
striction terms giving rise to nondeviatoric spontaneous strain in addition to the deviatoric 
one. For the order of magnitude estimate done below this difference is, however, irrelevant. 

LiNbOs possesses the symmetry i?3c in the high-temperature paraelectric phase and 
exhibits a 2nd order transition into the ferroelectric i?3c phase at 1460 K [82|, the transition 
being described by a one-component order parameter. 

Material parameters of BaTiOs, PbTiOs and LiNbOs are summarized in the Table 1. 

Here T^ is the melting point. In the case of BaTiOs and PbTiOs Tq is the temperature of 
the hrst order transition, while for LiNbOs this value yields the Curie point. Note tahe the 
scatter of the Kjc values for BaTiOs (Table 1 ) originates from experimental results obtained 
on ceramics with different properties (such as porosity, grain size, etc.). It gives rise to the 
corresponding spread of the estimates for BaTiOs parameters summarized in Table 2 . In the 
case of BaTiOs the last row yields the values of the phase transition line slopes corresponding 
to each its transition; II-III, III-IV and IV-V, 

The estimates following form the above material constants are collected in Table 2. 

In BaTiOs the slopes of all phase transition lines are negative. For this reason the 
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Table 2: Estimates of t 

le LPT characteristics 


AT* (K) 

at*/t^ 

Lf (nm) 

Vc (xlO® cm / s) 

BaTiOa II-III 

10=^ to 10^ 

1 to 10 

0.1 to 1 

1 

BaTiOa III-IV 

10^ to 10^ 

0.1 to 1 

1 

0.1 to 1 

BaTiOa IV-V 

10^ to 10^ 

0.1 to 1 

1 

0.1 to 1 

PbTiOs 

10^ 

1 

1 

1 

LiNbOa 

30 

0.01 

10 

- 


transformation zones only show up above the corresponding transition lines. Thus, during 
fracture of the cubic BaTiOa one will hnd at the tip a zone containing the tetragonal phase 
III embedded into the cubic matrix II. The above estimates show that this zone exists up to 
the temperatures of at least ~ 100 to ~ 1000 K above the temperature T//./// = 394 K of 
the bulk phase transition. 

The zone containing the orthorhombic phase IV embedded into the bulk tetragonal phase 
III should be found at the crack tip taking place above the bulk transition temperature 
Tin-IV- It should be especially noted that the temperature interval of existence of the 
orthorhombic zone, AT^ni-iv ~ 100 to 1000K, is greater than the ’’distance” between 
the transitions II-III and III-IV: Tn-iii — Tni-iv = HOK. Thus, the orthorhombic zone 
shows up in the whole region of existence of the tetragonal phase and will is detectable 
within the bulk cubic phase together with the tetragonal zone. In other words, a two-phase 
transformation zone should take place over the temperature Tn-iii- 

It is worth noting that the room temperature belongs to the temperature interval in 
which the tetragonal bulk phase III exists. It makes observation of the local phase transition 
phenomenon in the phase III of BaTiOa convenient. 

Finally, the zone with the ferroelectric, rhomboedral phase V should take place on the 
background of the ferroelectric, orthorhombic, bulk phase IV up to ~ 100 to 1000 K above 
the bulk transition temperature Tiv-v- One hnds Tiv—Ty = 84 K and ATnv-v ^ Tiv—Ty- 
This implies that the rhombohedral zone should be detectable within the whole domain of 
existence of the bulk, orthorhombic phase IV and as well as in the tetragonal matrix III. 
A proper description of such multiphase zones is, however, complex; we leave it for another 
paper. 

Since in PbTiOa one hnds fc < 0, it only exhibits the transformation zone during fracture 
of the cubic phase (i.e. above the transition line of the phase diagram). From its phase 
diagram m one can see that at the atmospheric pressure the cubic phase only exists at the 
temperature over about Tq ~ 850 K. If the PbTiOs fracture takes at temperatures between 
To and about 1000 K above Tq, our estimates predict that it will be followed by formation of 
the tetragonal transformation zone at the crack tip embedded into the matrix of the cubic 
phase. 

It should be noted that both in the case of BaTiOa and PbTiOa the ratio of the temper¬ 
ature shift, AT*, to the melting temperature, T^, is between ~ 0.1 and ~ 1. Since the whole 
phase diagram spans between OK and T^, this shows that the region of the zone existence 
covers a considerable part of or even the whole phase diagram. 

In principle, the slope, A;, of the phase diagram line may have any value, including a very 
large or a small one. The latter is the case of LiNbOa where it gives rise to a relatively small 
shift of ~ 10 K. 

5.3.2 Typical values of the temperature shift AT* 

It should be noted that the huge values of AT* ~ 100 to 1000 K are not only inherent for 
BaTiOs, PbTiOa. To argue that let us note the expression ([2]) for the temperature shift at 
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y = 0 can be written as: 


AT, ~ rJ/\kKjc)^/^ (58) 

where r^o = {g/ay^"^ is the order parameter correlation radius at the ’’distance” of 1 K from 
the transition line. The latter represents one of the most accessible parameters, since it can 
be extracted from the width of the X-ray spectrum peaks EZI typically exhibiting the value 
of Tco ~ 1 to 10 nmK^^^, as well as that there is a typical value of the slope of the phase 
diagram line fc ~ 1 to 10 K /kbar ~ (0.1 1) x 10“® Kcm^erg“^ [88], and the typical value 

of the fracture toughness of the inorganic solids is Kjc ~ 1 MPa m^/^ ~ 10® erg cm“®/^ [50] 
one hnds the typical temperature shift: 

AT, ~ 10^ to 10^ K 

Since the typical values of the melting point of inorganic solids is Tm ~ 10® K, one concludes 
that AT, typically covers a considerable part, if not the whole phase diagram above or below 
the line of the bulk phase transition. 

5.4 On difficulties and possibilities to detect a transformation pro¬ 
cess zone outside of the hysteresis region 

In the introduction we listed a number of materials for which LPT observation has been 
reported. A relatively small number of such materials others than those of the martensite 
type seems to contradict our main hndings. This is, however, only an apparent contradiction. 
Detailed inspection of the works cited above shows that in most of them the local phase 
transition has been detected by the analysis of the fracture surface available after the sample 
has been broken, the so-called, ’’post mortem” examination. In the case of the zirconia, for 
example, the fracture surface exhibited a layer of the monoclinic daughter phase about 1 pm 
thick on top of the tetragonal mother phase surviving a considerable time after fracture. This 
requires the zirconia to be deep within the hysteresis region of its phase diagram. Indeed, 
martensitic transformations in both martensite-austenite metals and zirconia exhibit wide 
hysteresis regions. This is, however, rare, for most solids the hysteresis does not exceed 
~ 10 K, while the second order transitions have no hysteresis at all. 

In contrast to the local phase transition within the hysteresis, the transformation zone 
outside of the hysteresis region is only present under stress. As soon as the stress is removed 
it immediately disappears. It cannot, therefore, be detected by the ’’post mortem” inspection 
of the fracture surface. 

Further, at a high temperature one may observe no zone at Kj = Kjc, but since 
AT, ~ the zone may show up at higher stress intensity factor, that is, at the tip 

of a propagating, rather than motionless crack. 

Detection of a local phase transition outside the hysteresis region is, therefore, a chal¬ 
lenging experimental task. 

5.5 The crack tip zone concept 

The notion ’’process zone” is already in use since long time. Its emergence reflects the 
understanding that the small domain in the immediate vicinity of the tip of a brittle crack 
should have special properties due to high stresses present there. Behavior of propagating 
cracks exhibits pronounced deviations from predictions of linear fracture mechanics [38] . It 
is generally believed that they can only be explained by mechanisms located within the 
process zone [10], [H]. Nevertheless, the content of this notion stays so far on an intuitive 
level. Below we outline the process zone concept from the physical point of view. 
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Equation of motion of the crack tip can be regarded as that of balance of forces, the 
driving force, K]/E, being balanced by (i) K]fj/E representing a kind of ’’dry” friction force 
and (ii) by additional resistance force generated within the process zone. Within the force 
balance concept the zone is the source of resistance of the solid to the crack propagation, 
since it is here that this additional force is generated. 

The concept of the process zone itself implies that it is possible to distinguish the solid 
inside from that outside of the zone using at least one physical property. Though in principle 
one can imagine a zone with a smooth, gradual variation of all its properties, we believe that 
the situation in which at least one physical property of the solid abruptly varies across the 
zone boundary is much more realistic and, hence, often met. 

We parametrized the differences of the solid properties inside the zone from those outside 
by the held, r]. Without the loss of generality we assume rj = 0 inside the zone, while 
vanishing outside. In this respect our point of view is akin to the popular phase held 
approach [HI]. In principle, this held may describe either an abrupt quantitative variation 
of some solid property, or its qualitative change. In the latter case the situation is usually 
qualihed as a phase transition, while rj is referred to as the order parameter. In the present 
paper we focus on this latter case. 

At present the terms ’’phase transition” and ’’order parameter” unify a crystal struc¬ 
ture variation of solids with bifurcations in non-linear systems, both equilibrium and non- 
equilibrium, such as e.g., bifurcations taking place during chemical reactions [H2]- It is tightly 
related to the fact that each of these transitions can be described by its inherent order pa¬ 
rameter (s). They exhibit a few classes of universalities generating the corresponding types 
of kinetic equations imposed on the order parameter(s), therefore, giving rise to different 
conhgurational forces acting on the zone boundary and, hence, exerted on the crack tip. 
The most often met universality class is associated with the order parameter obeying the 
Ginzburg-Landau-Khalatnikov equation [92], since it is related to the simplest bifurcation. 
This one is addressed in the present paper. The zone at the crack tip may, however, be 
associated with any of such phase transitions or their combinations. 

The long-wave elastic field at the tip of a brittle crack exhibits its own universal behavior 
expressed by the well-known small-scale approximation for the stress or strain field eu ~ 
^- 1/2 fjpjg latter universality combines with one of those mentioned above. 

We illustrate the above ideas within the example of classical structural phase transitions 
related to changes in lattice structures in crystals in response to the temperature and/or pres¬ 
sure variations. Structural phase transitions can be classihed according to their own classes 
of universalities. These are related hrst of all to the symmetry change taking place during 
the transition(s) [SS], [59], [HO]. The latter is manifested in (i) the number of components 
of the order parameter and (ii) structure of the Ginzburg-Landau-Khalatnikov equation. It, 
in particular, defines the form of interaction of the order parameter with other degrees of 
freedom including elastic ones. 

One observes in addition that the ranges of the values of material parameters in use 
are quite narrow. This allows one to determine typical numerical values of the parameters 
derived in our model. 

To conclude, within our approach the process zone is regarded as a domain which is 
unambiguously different from the bulk of the solid. The difference may either be quantitative 
(such as a prominent variation of at least one physical parameter), or qualitative (e.g. the 
difference in its symmetry, crystal structure or chemical composition). Such a variation 
(or their combination) takes place due to the high stress in the vicinity of the crack tip, 
as a consequence of the nonlinearity. Our approach is, further, based on the accounting 
for the universalities inherent both to the fracture mechanics and bifurcations of non-linear 
systems offering a way to classify zones according to the universality classes as well as on 
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the observation of their typical features. The zones related to the elastic nonlinearity [13], 
formation of the secondary cracks 001 , crack tip chemical reactions [93], electronic structure 
variation [93] as well as the transformation zone reported in the present paper £t into this 
general scheme. 


6 Appendix A: Exclusion of the acoustic variables 


Below we exclude acoustic degrees of freedom from the equation of motion as it has been 
proposed in the paper [IS]. Making use of f[T91 [T9]) one can express the displacement vector, 
Ui as 

Ui = uf\r) + A j (59) 

where is the displacement held created by the ’’undressed” crack (that is, the crack 

without any LPT) and Gij{r) is the Green function of the elastic solid with a cut. To the best 
of our knowledge the explicit form of such a Green function is unknown. We approximate 
it with the Green function of the inhnite elastically-isotropic body [66]. Passing to the 
reciprocal space under the integral in the representation 0591) and making use of the identity 


Gik{q)qjqk 


1 - 2(7 qiQj 
2p(l-a) q2 


(where Gjfc(q) is the Fourier-transform of the Green function) one hnds the strain held, 
eik{r) in the following form: 


Gfc(r) 



kl(l -2cr) 
2/i(l - a) 


^Q(q) exp(iqr) 


d?q 

{2nf 


(60) 


where £i°^(r) is the strain of the ’’undressed” crack and 


Q(q) = 2Ti5{q^) 


y) exp[i(gix -|- q 2 y)dxdy 


is the Fourier-image of ? 7 ^(r). The result fl60l) leads one to the following expression for the 
trace of the strain held: 

G*(r) = (r) + 2^1^ - 1 ) 

The hrst term in the right-hand sides of any of the expressions fl59|) . fl6n|) and fl6T|) describes 
the strain held generated by the ’’undressed” crack, while the second one yields the LPT 
contribution. Substitution of the strain trace dSnilHID into the free energy ([5l[6l[71[8|) yields 
the ehective free energy flTTI) with the factor /3 flT6|) in front of instead of /3q. 


7 Appendix B: Branching equation in the low-temperature 
phase 

In the low-temperature phase the branching equation is mostly convenient to obtain starting 
from the ehective free energy. The latter has the form: 

Fes = Fq + I 2 (q;*2 — «) ^2 + 3-^3^2 —^^2 (62) 
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It differs from the expression by the existence of the cubic term, where 






rdrdd 


( 63 ) 


Passing to dimensionless variables: r = P/R 21 V = vV^ one hnds: 


h{i^ 2 )= A{p,9)B{p,9)pdpd9 


(64) 


where 


A{p, 9) = exp {3 X 2^/^y'pcos(6*/2) — 3 x 2‘^^^p [1 + z/cos(6*)]} 




2 X 2 i/» - " 2 ) v^ - 2 X 2"''='cos(«/2)]"''^ 


The integrals I 2 and 1 4 are defined as in flT2]) . 

The integrand of I 3 becomes complex as soon as the expression (1 — i/|) y/p — 2 x 
2^/^ cos(6*/2) under the radical becomes negative. This takes place along the line 


Po{0) 


4 X 2^/^cos^(6'/2) 


(65) 


where ? 7 o(r) iHTIl turns into zero. At smaller values of r one finds p = 0. For this reason in 
A, 2,3 one should only integrate over p from po{9) to infinity, while the integration over 9 runs 
from —TT to TT. 

The integration has been done numerically using a standard NIntegrate routine of Mathe- 
matica 10.1 [70] employing an even-odd subdivision method with the local adaptive strategy. 
Below only the ratios / 2 / 4 //I and R/I 4 are used. These ratios and their fitting by simple 
functions: 

si(i/) = hh/lj ~ 0.59 + 3.131/2 - 9.50z/| + 9.00//^ - 3.23iyl ( 66 ) 

0 31 

S 2 M = h/h ~ 1-45 — 1.85z/2 -|- 2.44i/| -|-^ (67) 

^ ^ ^ l-0.98i/| ^ ^ 

The effective free energy fl62]) has a cubic term. Since it is positive, one finds that the 
left minimum of the free energy fl6^ is more pronounced. Analogously to the static state 
one concludes that the solution of the branching equation should be chosen that corresponds 
to this deeper minimum, that is, the negative one. Let us note that in a general case the 
sign of 6 p is opposite to the one of Pq. Should we have chosen a negative po, we will get the 
positive sign for 6 p. 

This solution of the branching equation takes the form: 


6 


S2{u2)B^A 

gl/6f]l/2 



1 + Si{u2) 


1 + 2 


glA 

54/3 




a — 


( 68 ) 


8 Appendix C. Simulation: technical details 

To perform simulations we used the software COMSOL 4b. Equations have been simulated 
in a half-plane t/ > 0. A semi-circular domain has been defined with the diameter, D = 50. 
By trial and error we find that it is large enough, to let the solution vanish well far from 
the domain boundary. The initial mesh size of 5 has been chosen, but the adaptive mesh 
refinement option has been further used to automatically refine the mesh as appropriate. 
The no-flux boundary condition has been set at the boundary y = 0 and the condition u = 0 
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Pseudo-time, tp^ 

Figure 12: Illustration of the converging of the calculations with the psendo-time showing 
the convergence of the amplitudes of the rescaled order parameter, Umax, at different q values. 
Filled disks: q = 0.3, Filed squares: 0.32, Filed diamonds: 0.34, vertex-np triangles: 0.36, 
vertex-down triangles: 0.38, open circles: 0.39, empty sqnares: 0.41. Note that in the case 
of g = 0.39 (point-down triangles) the Mmax(t) dependence still exhibits a slope and a for the 
satisfactory convergence a longer process was nsed (not shown). 


at the rest of its bonndary. A straightforward simnlation of the static equation fl5^ with 
such boundary conditions, however, only retnrns the trivial solntion u = 0 at any value 
of the control parameter q. To avoid this instead of flS5]l we introduced a pseudo-dynamic 
equation: 


du ^ .in du 

-= Aiu + 2^/V— 

UX\ 


dt 


ps 


{xl + + xi 


1 / 2 ' 


{xl + yl + ef^"^ 


u 


u 


( 69 ) 


1 1/2 

where u = u(a;i, gi, tps), fps is the pseudo-time and {xl -|- + xi / {xf + yf + is 


equal to cos(6'/2)/r}'^^, regularized in the vicinity of the point ri = 0 by a small parameter 
£ = 0.0001. Stable solutions of the static equation fl53|l represent Fxed points of the dynamic 
system fl69|l . As initial condition we used a smoothed step function, only unequal to zero in 
a vicinity of the point (0, 0). 

The dynamic system has been solved using the direct MUMPS solver with the BDF 
time stepping. The convergence of the solution to its Fxed point has been controlled by 
the behavior of the Umax, the maximum value of the function u{xi^yi,tps). By trials we 
found that 700 pseudo-time steps ensure a good convergence, though sometimes it has been 
necessary to keep the process as long as 3000 steps. Figure [12] shows the example of such a 
convergence study for a number of simulations in which all parameters except q were Fxed, 
while q varied. 

One can see that far from the bifurcation point the convergence takes place well before 
700 pseudo-time steps are done. As it can be expected, the situation is diFerent in the close 
vicinity of the bifurcation (g = 0.38 and 0.39 corresponding to the vertex-down triangles 
and open circles in Fig. [T2]) . Even here 700 pseudo-time steps guarantee a rather reliable 
convergence. 
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